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Spherical harmonics analysis of Fermi gamma-ray data 
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We argue that the decomposition of gamma-ray maps in spherical harmonics is a sensitive tool to 
study dark matter (DM) annihilation or decay in the main Galactic halo of the Milky Way. Using 
the spherical harmonic decomposition in a window excluding the Galactic plane, we show for one 
year of Fermi data that adding a spherical template (such as a line-of-sight DM annihilation profile) 
to an astrophysical background significantly reduces x^ oi the fit to the data. In some energy bins 
the significance of this DM-like fraction is above three sigma. This can be viewed as a hint of DM 
annihilation signal, although astrophysical sources cannot be ruled out at this moment. We use 
the derived DM fraction as a conservative upper limit on DM annihilation signal. In the case of 
bb annihilation channel the limits are about a factor of two less constraining than the limits from 
dwarf galaxies. The uncertainty of our method is dominated by systematics related to modeling 
the astrophysical background. We show that with one year of Fermi data the statistical sensitivity 
would be sufficient to detect DM annihilation with thermal freeze out cross section for masses below 
100 GeV. 

PACS numbers: 95.85.Pw, 95.55.Ka, 98.70.Rz, 95.35.+d 
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MOTIVATION 



Out of all indirect searches for dark matter 
gamma-rays are probably the most "direct" 
Charged particles, such as positrons and antiprotons, 
are deflected in the Galactic magnetic field. The infor- 
mation about their source is lost and only anomalies in 
the spectrum may signal the presence of DM. Most of 
gamma-rays propagate freely inside the Galaxy and, 
together with the spectrum, they carry information 
about the morphology of the source. This property 
may be crucial in separating a DM signal from astro- 
physical backgrounds, e.g., ^3,]. 
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From cosmological simulations y42| , we expect that 
cold dark matter in our Galaxy has formed a nearly 
spherical halo with density growing toward the Galac- 
tic center (GC). Thus, DM annihilation or decay may 
be a source of gamma-rays with a spherical shape 
peaked at the GC, in addition to astrophysical and 
extra-galactic sources. 

In this paper, we study the contribution from the 
main spherical halo ignoring DM substructure. In or- 
der to minimize the astrophysical fiux, we mask the 
Galactic plane and resolved gamma-ray point sources. 
The problem is that at high latitudes a possible DM 
annihilation signal is relatively smooth and most prob- 
ably subdominant to Galactic and extragalactic diffuse 
emission. In the paper we propose to use the spher- 
ical harmonics decomposition of gamma-ray data to 
search for DM annihilation or decay. The contribution 
of a smooth signal with small amplitude is maximal 
for spherical harmonics with small angular numbers i. 
Consequently, the Galactic DM signal away from the 
GC may contribute most significantly to small £ har- 



monies, while its contribution to large I harmonics can 
be neglected compared to the Poisson noise. 

Spherical harmonics decomposition has several ad- 
vantages compared to template fitting in coordinate 
space: 

1. Organization of data: small i harmonics carry 
all the information about the large-scale distri- 
bution of sources, while large £ harmonics are 
dominated by the Poisson noise. Spherical har- 
monics decomposition is a linear transformation 
that has no information loss, but only relevant 
information for large-scale distributions is used 
in fitting. 

2. Universality: small i harmonics are insensitive 
to the resolution of pixel maps (for sufficiently 
small pixel sizes). In particular, the templates 
may have different resolution and do not need to 
be brought to the same pixel size as in the case 
of coordinate space fitting. The x^ is also inde- 
pendent of resolution, while in coordinate space 
the absolute likelihood depends on the number of 
pixels. 

3. Linearity and stability: both the transforma- 
tion of data from coordinate space and the fitting 
in spherical harmonics space are linear operations 
(the x^ has the usual quadratic form). Thus, a 
nonlinear Poisson likelihood in coordinate space 
is substituted by a combination of two linear op- 
erations in spherical harmonics space. This may 
be useful for stability of the fitting procedure in 
the case of small numbers of photons in pixels: 
the Poisson probability is undefined for nonpos- 
itive expected numbers of photons, while small 
negative expected numbers should not a be prob- 
lem: it simply means that the template is not 
perfect and we over subtract this template to fit 
data somewhere else. 

4. Symmetry: in some cases spherical harmonics 
may be useful to focus on a part of data with a 
particular symmetry in mind. In this paper we 
will use the spherical symmetry of dark matter 
distribution around the Galactic center. If we 
point the z-axis toward the GC, then dark matter 



contributes only to Yi„i harmonics with ?7i = 0, 
i.e., we may select only a/o modes in fitting. 

A computational algorithm for fitting in spherical 
harmonics space is straightforward but there are a few 
things to keep in mind. First, the Yim's are not orthog- 
onal in a window on the sphere. The corresponding 
spherical modes, aim's, are still independent but corre- 
lated. As a result, their covariance matrix is not diago- 
nal. In general, this may render the computations un- 
feasible, unless one uses some special techniques (such 
as the Gabor transform for the power spectrum [Sj). 
We will use only a few Yjo harmonics corresponding to 
the largest scales. The corresponding covariance ma- 
trix is relatively small and can be easily computed. 

The choice of the astrophysical background model is 
a more conceptual problem. A thorough solution of 
this problem can be quite complicated and we will not 
discuss it here. The main purpose of our work will 
be to illustrate the method of the spherical harmonics 
transform in the analysis of gamma-rays. As a toy 
model for the astrophysical background we will use the 
gamma-ray distribution in a low-energy bin, since we 
expect that the DM contribution to the spectrum is 
insignificant at low energies. 

The paper is organized as follows. In Section HIl we 
describe an algorithm of fitting templates in spherical 
harmonics space. In Section IlIII we apply this method 
for Fermi gamma-ray data. We compare two cases. 
In the first case we use two templates: a low-energy 
bin and an isotropic distribution. In the second case 
we also add a distribution of photons with a spherical 
symmetry around the GC. We find that the residual for 
the three-template models has a much better x^ than 
the residual for the two-template model. In Section 
IIVI we find the best-fit energy spectrum of the fluxes 
assuming a power-law dependence on energy. We also 
put constraints on DM annihilation in bb. Section |V] 
has conclusions. 

There are three appendixes. In Appendix [A1 we cal- 
culate the covariance matrix for spherical harmonics 
defined on a window in the sphere. In Appendix IbI we 
check the fitting algorithm with a Monte Carlo simu- 
lation. In Appendix [C] we discuss the contribution to 
the angular power spectrum from point sources. 



II. 



METHOD 



In this section we describe a general method of tem- 
plate fitting in spherical harmonics space. In the next 
section we apply this method to the Fermi gamma-ray 
data to search for DM annihilation in the Milky Way 
halo. 

In general, an algorithm will contain several steps: 

1. Choose a mask (for instance, one can mask the 
Galactic plane and point sources). 

2. Find the spherical harmonics decomposition of 
the data outside of the masked region, aim- 

3. Find the covariance matrix for the spherical har- 
monics coefficients, Cov{aim,ai'm')- 

4. Formulate a model for the gamma-ray distri- 
bution as a function of parameters a and find 
the corresponding decomposition into spherical 
modes bim (a) outside of the mask. 

5. Find the best-fit model parameters a* by min- 
imizing a x^j where we use the full covariance 
matrix instead of a^ due to a nontrivial correla- 
tion of the spherical modes on a window in the 
sphere (Equation (jASp V 

In the remaining part we will mostly introduce the 
notations that will be necessary in interpreting the re- 
sults of data analysis in the next section. The mathe- 
matical details can be found in Appendix [A] 

In the calculations of spherical harmonics decompo- 
sition, it is convenient to use a pixelation of the sphere 
(we use HEALPix [9|). We will consider some energy 
bins Ei and denote by np{Ei) the number of photons 
inside the energy bin i in a pixel p — 1, . . . , A^pix- For 
clarity, in the following formulas we will suppress the 
dependence on energy. 

Define the photon density as p(7p) = ^p, where 
Srip is the size and jp is the center of the pixel p. We 
put Up = p^jp) — 0, if the center of the pixel is inside 
the mask. The spherical harmonics transform of the 
gamma-ray data is 

aim = / Y;^h)p{i)dn « Y. YLiipW (1) 



In Appendix |A] we show that the covariance matrix has 
the following simple expression 



Gov {aim, ai'm') 
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Yirm'ilp)Ylmilp)np. (2) 



Let us now describe the parametrization of the mod- 
els. In general, the shape of the model fluxes can de- 
pend on some parameters. In this paper we will focus 
on the template fitting, where the shape of the tem- 
plates for various sources is fixed and the only variable 
parameters are the normalizations in every energy bin. 

For each template a in every energy bin, the variable 
parameter will be the average flux inside the window 
F" in units of ■ 



„ ,, ■) . In order to find the number 

CjC V cin^ s sr 

of photons in a pixel p from a template a, one needs 
to multiply F" by the Fermi exposure Ep times the 
probability distribution function (PDF) p" that carries 
the information about the shape of the template a, 
times the size of the window ftw, times the size of 
the energy bin AE. The number of photons from the 
template a in a pixel p is 



1^ = i^" • {Ep ■pf.-nw AE), 



(3) 



where the PDF is normalized as J^pPp — 1- Let us 
denote the spherical harmonics decomposition of the 
function in the parenthesis as 

vL = E ^™ (>) (^P • Pp • f^w AE) , (4) 

p=i 

then the spherical harmonics decomposition of a linear 
sum of sources is 



blmiFn=J2P' 



•an 



(5) 



The best-fit fluxes F^ can be found from minimizing 
the x^ in Equation (jASp . 

The uncertainty of the model parameters F" around 
the minimum of x" can be estimated from the Hessian 
matrix 



Ha[3 = 



dF"dFl^ 



(6) 



F"=F° 



In particular, the variance of the model fluxes can be 
estimated as 



Var(F") = H-\ 



(7) 



III. DATA ANALYSIS 



2. Rotate the z-axis to point toward the GC. 



A. Data selection 

We consider 13 months of Fermi gamma-ray data 
(August 4, 2008 to August 25, 2009) that belong to the 
"diffuse class" (Class 3) of the LAT pipeline. We ex- 
clude the data beyond zenith angles of 105° due to sig- 
nificant contamination from atmospheric gamma-rays. 
We also exclude the data taken over the South Atlantic 
Anomaly (SAA) and mask the point sources detected 
by Fermi [10[. 

Most of the time we will use the gamma-rays between 
1 GeV and 300 GeV which we separate in 10 exponen- 
tial energy bins between 1 GeV and 100 GeV plus an 
extra energy bin between 100 GeV and 300 GeV. We 
mask the pixels centered within 10° from the Galactic 
plane. We also mask all pixels that either contain a 
gamma-ray point source or if the boundary of the pixel 
is within 68% containment angle Ki 0.7° at E ^ 1 GeV 
11| to a gamma-ray point source. 

The interpretation of spherical harmonics decompo- 
sition of a DM model is simplest in the coordinate sys- 
tem where the z-axis points toward the GC (at odds 
with the standard Galactic coordinates in which the z- 
axis points toward the Galactic North pole) . We choose 
the a;-axis pointing toward the Galactic South pole. If 
we were considering all data without masking, then DM 
would contribute only to Yiq harmonics due to a rota- 
tional symmetry round the new z-axis. In the presence 
of a mask, DM contributes to all spherical harmonics, 
but its contribution to Yiq modes is still maximal and 
we will restrict our attention to these modes for sim- 
plicity of analysis. 

Pixelation of the data and spherical harmonics de- 
composition is performed with healpy, the python ver- 
sion of HE ALPix Q[35|. 

An example of a gamma-ray counts map for the en- 
ergy bin between 1 GeV and 1.6 GeV with the z-axis 
pointing toward the GC can be found in Figure [T] 

Summary of data selection and model parameters: 

1. Mask the gamma-ray point sources and the 
Galactic plane within |6| < 10°. 



3. Consider Yim harmonics with ^ < 15 and m — 0. 

4. We choose HEALPix parameter nside = 32 (cor- 
responding to pixel size of about 2°). 

5. Astrophysics template energy bin: 1 GeV < E < 
1.6 GeV. 

These are the data selection parameters for the "main" 
example that we consider in Section IIIIBI In Section 
nil CI we consider the effect of varying these parameters. 



B. Bin to bin fitting 

We use the energy bin between 1 GeV and 1.6 GeV as 
a template for the Galactic astrophysical emission as- 
suming that the gamma-ray emission at these energies 
is dominated by the Galactic cosmic ray production 
through the 7r° decay. One of the possible limitations of 
this template is the inverse Compton scattering (ICS) 
component of gamma-rays. Relative contribution from 
the ICS photons increases with energy and may become 
comparable to tt'^ photons above 10 GeV [ij]. As a re- 
sult, a low-energy bin template may underestimate the 
ICS component at higher energies. Currently, there is 
no generally accepted model for the ICS emission (see, 
e.g., the caveats section in the description of the Fermi 
diffuse background Imodell [sef ) . We will treat the ICS 
component as a systematic uncertainty in the current 
analysis. 

We consider an isotropic template as a model for 
the extragalactic emission. We will also consider sev- 
eral templates with a spherical distribution around the 
GC. Our main working example will be the line-of-sight 
DM annihilation in Navarro, Frenk, and White (NFW) 
profile 0. The NFW profile is 



PDMir) oc 



1 



r, (l + r/r,)2 



(8) 



where the scale parameter Ts = 20 kpc 13|, [ij] . The 
window angle |6| > 10° corresponds to distances r > 
1.5 kpc from the GC. At these distances the NFW pro- 
file is similar to less cuspy profiles, such as the Einasto 



Masked map, 1.0 - 1.6 GeV 




FIG. 1: Counts in pixels for 1 - 1.6 GeV energy bin. The Galactic center is in the north pole, the anti-center is in the south 
pole (southern hemisphere is in the center). We mask \b\ < 10° and gamma-ray point sources detected by Fermi [10|. 
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FIG. 2: An example of the aim fitting procedure for m = and £ = 0, . . . , 15 in the 6.3 GeV to 10 GeV energy bin. The 
harmonic decomposition in a window is defined in Section [IT] (there is an additional HEALPix normalization factor 47r/A'pix 
with respect to Equation H}). The low-energy bin aim's are given by the gamma-ray data in 1 GeV to 1.6 GeV energy 
bin. The isotropic template aim's are nonzero for ^ > due to the spherical harmonics decomposition in a window. Large 
fluctuations of data aim's are mostly due to symmetry of the mask and the fact that the Galactic flux is stronger near the 
Galactic plane (notice, that the modes with even i are significantly larger than the modes with odd ^'s). The noise level is 
given by the square root of the diagonal elements in the covariance matrix in Equation ^ . The x^ for the two-template 
fit on the left is 80 for 14 dof, while the three-template fit on the right has x^ = 16 for 13 dof. The x^ in other energy bins 
can be found in Section Till D I where we also compare the NEW annihilation profile with other profiles. 



profile [15|. The annihilation signal is proportional to 
Pdmj '^f-i Equation [T51 

As a null hypothesis we take a combination of two 
templates: an astrophysical template modeled by a 



low-energy bin and an isotropic flux. We compare this 
model with the three-template model, where we also 
add a template with spherical symmetry around the 
GC. In fitting, we decompose the templates and the 
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FIG. 3: Left plot: two-template model fitting. Right plot: three-template model fitting. Points represent fitting in every 
energy bin independently. Lines represent best-fit with power-law spectra discussed later in Section ITVl The first bin is 
"singular" since we use the data in this bin as one of the templates. The points associated with the NFW annihilation 
template are significantly above zero, i.e., there is a significant signal that grows toward the GC, but, as we show in 
Figure|6l other profiles with rotational symmetry around the GC give similar improvement of x^, i.e., current data are not 
sufficient to distinguish different profiles and understand the nature of the signal. In Section IIVBI we will use the points 
associated with the NFW annihilation profile to put conservative upper limits on annihilation cross-section. 



data in every energy bin into spherical harmonics ac- 
cording to the algorithm in Section [III The fluxes cor- 
responding to the templates are found by minimizing 
the x^ in Equation (jA5| . An example of fitting the 
data in an energy bin between 6.3 and 10 GeV using 
the spherical harmonics is presented in Figure [H The 
results of fitting in every energy bin are shown in Fig- 
ure [H 

The best-fit value of the flux corresponding to the 
spherical template can be used to put a conserva- 
tive upper bound on DM annihilation into a pair of 
monochromatic photons (the DM line signal) . The up- 
per bounds obtained this way are an order of magni- 
tude less constraining than the limits on monochro- 
matic gam ma-rays obtained by the Fermi Collabo- 



ration [16|. The spherical harmonics decomposition 
method has a better performance for signals with 
smooth energy dependence. In Section IIVBI we find 
conservative Umits on bb DM annihilation only a factor 
of a few less constraining than the limits from dwarf 
spheroidal galaxies [17.] . 



C. Variation of parameters 

In this subsection we check the robustness of spheri- 
cal harmonics decomposition with respect to variation 
of data selection and model parameters. The main 
model is given by the flux associated with NFW anni- 
hilation template found in Section IIIIBI with parame- 
ters defined in Section UlIAl We consider the following 
variations of parameters (Figure 2]): 

1. Different energy bins for the astrophysics tem- 
plate: 0.7 GeV <E <l GeV and 2.5 GeV <E < 
4 GeV. The contribution above 10 GcV is consis- 
tent with the main model. 

2. The DM flux does not change significantly if we 
shrink the window size to |6| < 5°, change the 
number of harmonics to i'max = 10, or change 
the resolution to nside = 64 (pixel size of about 

V). 

3. Separating the northern and the southern hcmi- 
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FIG. 4: Dependence of flux associated with the NFW an- 
nihilation template on model parameters. The main model 
and the variations are described in Sections lIIIBI and llllCI 



spheres: the contribution of the spherical tem- 
plate in the north is significantly smaller than in 
the south. This may be due to a stronger as- 
trophysical flux in the north (compare with the 
Schlegel, Finkbeiner, Davis (SFD) dust template 



18| that is believed to trace the tt'^ production of 



gamma-rays 



H 



We performed other variations, such as increasing the 
window size and increasing ^max (not shown in Figure 
131). We find that the variation of parameters does not 
change the results significantly. 



D. Different spherical profiles 

In this section we compare different profiles with 
spherical symmetry. Together with the NFW DM an- 
nihilation profile ex Pnfw studied in Section [ill Bl we 
consider NFW DM decay c>c pnfw, a profile oc I/r^ 
corresponding to the distribution of mass in the stellar 

(r is the distance from 
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halo of the Milky Way 

the GC), and a bivariate Gaussian profile with ui = 15° 

and (76 — 25° studied in [l9|. 
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FIG. 5: Best-fit values of the flux associated with different 
spherical proflles. The values are found by fltting a three- 
template model (low-energy bin, isotropic, and spherical 
templates) to the data. The case of NFW annihilation is 
discussed in detail in Section Till Bl 



The best-fit fiux values associated with these profiles 
are shown in Figure [5] In Figure |6] we compare the x^ 
for these model with the null hypothesis (a low-energy 
bin template plus an isotropic flux) from Section ITlIB I 
We find that any of the spherical profiles give a signifi- 
cant improvement of x^ compared to the two-template 
model. The stellar halo profile has smaller y^ than the 
other three profiles. A possible astrophysical source of 
high energy gamma-rays at high latitudes could be a 
population of millisecond pulsars 23|, 
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In Figure[7]we compare our fit of the bivariate Gaus- 
sian halo with the calculation in [19| . There is a general 
agreement above ~ 2 GeV. Our error bars are larger 
than the errors in 19|, possibly due to the usage of a 
subset of spherical modes with to = rather than all 
spherical harmonics. 



IV. 



ENERGY SPECTRUM 



In this section we fit the fiuxes found in the pre- 
vious section by some general energy spectra. We 
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FIG. 6: Goodness of fit of the models presented in Figure 
[5] compared to the astrophysical model. The number of 
degrees of freedom (DOF) for the three-template model is 
13 (sixteen aio modes minus three parameters associated 
with the normalization of the templates) . 



assume that there are three main contributions to 
the gamma-ray flux: Galactic astrophysics emission, 
isotropic emission (extragalactic plus a possible con- 
tamination from misidentified cosmic rays), and an ad- 
ditional spherically symmetric flux. We flt the Galactic 
and the isotropic fluxes by power-law spectra. For the 
spherical template, we use a power-law with an expo- 
nential cutoff and an energy spectrum for DM annihi- 
lating into bb. 



A. Power-law energy spectra 

In the previous section we have used a low-energy 
bin as a template for the Galactic astrophysical com- 
ponent. In reality the flux in the flrst bin Eq is a sum 
of all components. As a result there is a nontrivial re- 
lation between the fluxes associated to templates, and 
the intrinsic fluxes. Let us denote the intrinsic Galactic 
component of the flux by <I'g(_E), the intrinsic isotropic 
component by $i(£'), and a spherical ("dark matter") 
component by ^d{E). The flux associated to the low- 
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FIG. 7: Comparison with the gamma-ray haze spectrum in 
[131 . Our flux is the same as in Figure [5] but avera ged over 
the window I e (-15°, 15°) and h G (-30°, -10°) [lil. 



energy bin will be denoted as F^{E), the flux associated 
to isotropic template as Fi{E), and the flux associated 
to spherical templates as Fd{E). 
The flux in the first energy bin is 



F,{Eo) = <^g{Eo) + $i(£;o) + $d(-Bo) 



(9) 



i.e., the corresponding template has contributions from 
the actual Galactic photons, from the isotropic dis- 
tribution, and, possibly, from an additional spherical 
component. Consequently, a fraction of isotropic and 
DM photons is included in the flux corresponding to 
the "astrophysical" low-energy bin template at all en- 
ergy bins. 

In this subsection, we consider the following 
parametrization of intrinsic fluxes 

E 

e'o. 



^. = $.n N^ 



gO 



*i = $in N^ 



E 

Eo 



(10) 
(11) 



*d = *do I ^^ 
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E'o 



(12) 



If we assume that the Galactic gamma-rays provide the 
most signiflcant contribution to the astro template at 



Model 


Jig 


ni 


rid 


Ecut (GeV) 


Astrophysical 


2.62 ±0.01 


1.49 ±0.03 


- 


- 


Decay 


2.79 ±0.02 


1.68 ±0.04 


2.13 ±0.05 


142 


Annihilation 


2.83 ±0.02 


1.88 ±0.03 


2.01 ±0.05 


137 


Stellar halo 


2.87 ±0.02 


1.89 ±0.04 


2.05 ±0.05 


152 


Gaussian halo 


2.84 ±0.02 


1.86 ±0.04 


1.98 ±0.05 


184 



TABLE I: Fitting the fluxes associated with difi^erent templates by power-law intrinsic Galactic, isotropic, and spherically 
symmetric fluxes. The parametrization is given in Equations (|13p . (I14p . and (|15p . The profiles are defined in Section [11101 



Eq, then we can expect that the astro template flux Fg^ 
has the same power-law index as the intrinsic Galac- 
tic flux $g and the following parametrization of i^a is 
reasonable: 



F,{E) = ($go + $io + $do) 



(13) 



The fluxes for the isotropic and DM templates are equal 
to the intrinsic fluxes minus the contribution to the 
astro template, 



FiiE) = $io 



$i, 



(14) 



F^{E) = $do f ;^j e ^^ - $do ( ;^ j (15) 

In order to find the parameters of the energy spec- 
tra, we use the fluxes found in Section IIIIBI in every 
energy bin as "data points" . The error bars are derived 
from Equations (HH), dH]), and (US]). The models of the 
fluxes are parameterized in Equations (fTO|). (fTT |) . (|12|) . 
The best-fit indices and the cutoff are presented in Ta- 
ble H 

The index rig « 2.8 for the Galactic component is 
consistent with the pion production of gamma-rays. 
The index n; « 1.9 for the isotropic flux is harder than 
the typical indices n — 2.2 — 2.5 for extragalactic diffuse 
background or Active Galactic Nuclei (AGN) spectra 
121 125| . This discrepancy is most probably due to an 
isotropic energy dependent contamination from cosmic 
rays (CR). A model for the CR background in diffuse 
class events can be found in Figure 1 of Rcf. [12:]. The 
corresponding spectrum is rather hard, ~ E~^, reach- 
ing '^ 0(1) fraction of the total flux around 100 GeV 
(compare Figure 1 and Table 1 in [12^]). 



Gamma-ray flux with an index n^ ~ 2 can be 
obtained by inverse Compton scattering of interstel- 
lar radiation photons and a population of electrons 
with a spectrum r^ E~^. A break at a few hundred 
GeV can be explained by transition between Thomp- 
son and Klein-Nishina scattering for star-light photons, 
£^brcak ^ [rrieC^Y / hv . If wc attribute the spherical sig- 
nal with the ICS photons, then the question would be 
to find a spherically symmetric source of high energy 
electrons at high latitudes. 

B. Limits on DM annihilation 

In this subsection we use the flux associated with the 
NFW annihilation template derived in Section UlIBI to 
put conservative upper limits on the rate of DM an- 
nihilation in the Milky Way halo. Assuming bb an- 
nihilation channel, we find the best-fit DM mass and 
annihilation cross section. In the analysis, we use the 
prompt gamma-rays emitted by the decay of bb, the 
corresponding spectrum is found with the help of the 
Pythia generator [26|, [27 1. The best- fit DM parame- 



ters are subject to significant systematic uncertainties 
due to modeling of astrophysical emission. The up- 
per limits, on the other hand, are rather robust. They 
only depend on the DM profile and on the annihilation 
channel, e.g., bb, W^W^ , etc. 

Let (av) denote the DM annihilation cross sec- 
tion. In this paper we will consider only the prompt 
gamma-ray emission from DM annihilation. The fiux 
of gamma-rays from DM annihilation per steradian at 
an angle 9 from the GC is 



Fbm{E,9) — 



1 



{av 



dN-y 



SvrAfgj^^ dE 



plM{r)dR (16) 
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FIG. 8: Left plot: upper bound on DM bb annihilation in an NFW profile in comparison with the dwarf spheroidal galaxies 
limit [17|. The bound is derived by using the NFW annihilation flux in Figure |3] as an upper limit. The statistical 
sensitivity is derived from the error bars for NFW flux in Figure [3] Right plot: best bb DM annihilation fit to NFW flux 
points in Figure [S] Contours represent Ax^. Limits from below are subject to systematic uncertainty due to modeling of 
astrophysical fluxes. Limits from above may depend slightly on DM proflle and annihilation channel. 



dN., 



where —j-^ is an average spectrum of prompt gamma- 
rays per annihilation event, R is the distance along 
the line-of-sight, and r is the distance from GC, r'^ — 



i?2 



Rfi + 2RRo cos ( 



Pdmo = 0.4 GeVcm 



We assume local DM density 

ii-Eo|. 



We parameterize the flux from DM annihilation by 
the DM mass and annihilation cross section. The result 
of fitting the corresponding energy density of the flux 
to the best-fit fluxes associated with the NFW annihila- 
tion template from Section llll Bl is shown in Figure|8]on 
the right. There are signiflcant systematic uncertain- 
ties, e.g., a distribution of inverse Compton photons, 
in proving the existence of a DM annihilation signal. 
One can nevertheless put upper limits on DM annihila- 
tion, provided that the flux from DM cannot be larger 
than the signal correlated with the NFW annihilation 
template. 

In Figure [5] on the left for every DM mass Mdm we 
find the best-flt annihilation cross section which gives 
the upper limit on DM annihilation. Statistical uncer- 
tainty in this case is an order of magnitude smaller than 
the limit itself. This uncertainty provides an estimate 



of statistical sensitivity of our method for the DM an- 
nihilation signal. For Mdm ^ 150 GeV this sensitivity 
is sufficient to detect DM annihilating with freeze-out 
cross section. 

In spherical harmonics analysis the uncertainty is 
dominated by systematics while statistical uncertainty 
is rather small. This makes spherical harmonics a com- 
plimentary tool to the searches of DM annihilation in 
dwarf spheroidal galaxies [iTJ, |3l| , where the systematic 
uncertainties are small and the limits are dominated by 
statistics. 



V. CONCLUSIONS 

In this paper we argue that the spherical harmon- 
ics decomposition is a convenient tool to study the 
large-scale distribution of gamma-rays such as a possi- 
ble contribution from DM annihilation or decay. The 
key points of this approach are the use of the spherical 
harmonics decomposition in a window that eliminates 
most of the known astrophysical sources and the choice 
of the coordinate system appropriate for the symme- 
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tries of DM distribution. In Appendix [B] we show that 
in a test with 10'' randomly generated photons, a 4% 
fraction of gamma-rays coming from DM annihilation 
can be detected with a five sigma significance. 

One of the main advantages of the spherical harmon- 
ics decomposition compared to the analysis in coordi- 
nate space is an efficient organization of the data. Con- 
sider, as an example, the top right plot in Figure 10. 
The harmonics with i! ^ 20 are dominated by the large- 
scale distribution of gamma-rays. The harmonics be- 
tween i ^ 20 and £ ^ 200 are dominated by the contri- 
bution from point sources, while the harmonics above 
£ '^ 200 are dominated by the noise (either physical 
noise due to insufficient number of photons or the Point 
Spread Function (PSF) of the instrument). Thus one 
can immediately separate the data that carry some in- 
formation about the large-scale distribution of gamma- 
rays from the harmonics dominated by the noise. 

Another approach in finding DM signatures in 
gamma-rays actively discussed in the literature is to 
search for featur es in the power spectrum due to DM 
subhalos 
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34| . Our method is complimentary to 



this, since we look for the signature of the main halo 
at small i, whereas DM subhalos usually contribute at 
£ > 100. We also believe that with the current vol- 
ume of data our approach is more advantageous since 
£ > 100 harmonics are dominated by the noise and 
much more data will be necessary to separate a sig- 
nificant signal, whereas for small £ the current data is 
enough to overcome the noise level for energies up to 
30 - 50 GeV. 

Our method already enables us to argue that there 
is a significant spherical distribution of photons in ad- 
dition to an astrophysical flux, which we model by tak- 
ing the data in a low-energy bin as a template, plus an 
isotropic distribution of photons. We compare several 
profiles with a spherical symmetry around the GC and 
find that below ~ 50 GeV "stellar" halo ex 1/r^ has a 
slightly better x^ than DM annihilation, DM decay, or 
a bivariate Gaussian distribution found bv |19|. 

We also use the flux associated with the NFW anni- 
hilation template to put upper limits on DM annihila- 
tion into bb. The derived limits are a factor of a few less 
stringent than the limits from dwarf spheroidal galax- 



ies. The uncertainty of our method is dominated by 
systematics while the dwarf spheroidal galaxies (dSph) 
method is dominated by statistics. One of the advan- 
tages of using the DM annihilation in the Milky Way 
halo versus the annihilation in dSph is the ability to 
use the Fermi data to put stronger constraints on DM 
annihilation, even after the Fermi LAT stops collecting 
data, by reducing the systematic uncertainty related to 
modeling the astrophysics backgrounds. 
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Appendix A: Spherical harmonics covariance 

In this appendix we provide details on the calculation 
of covariance matrix for spherical harmonics defined in 
a window on the sphere. 

Consider a pixelation of the sphere. We will assume 
that the pixel size is sufficiently small compared to the 
angular size of interest, but it is sufficiently large so 
that the photons in different pixels are uncorrelated. 
Denote by Up the number of photons in a pixel p, let 
7p be the center of p, and 6flp be the area of p. We 
can define a discrete photon density function pp = jjy- . 
The spherical harmonics transform of the density func- 
tion is 

aim = / Y:^i^)p{^)dn « Yl y*mhp>P^ (Al) 

•' p=l 

The spherical functions are defined as 



Y, 



Im \ 



"'-i/^i^''"'"""''""' 



(A2) 



12 



The spherical harmonics are orthogonal on the sphere 
but on a window in the sphere they are not. As a result, 
we expect correlation between different harmonics. 
The covariance matrix of a;,„ 's is 

ClmJ'm' = (aL«/'m') ~ ('^rm)(«/'m') (^3) 

p 

where we have used that the numbers of photons at 
different points are not correlated {iipUpi) — {np){npi). 
For a Poisson distribution (n^) — (rip)^ — (rip). In 
a particular realization of the photon map, our best 
estimate of (up) is the actual number of photons in 
this pixel Up. Consequently, the covariance matrix can 
be estimated as 



100 random tests, input DM fraction = 0.04 



t^im, I'' 



P 



Yl*m'ilp)Ylr>ihp)np- 



(A4) 



Spherical harmonics coefficients a;,„ together with their 
covariance matrix provide the data necessary to for- 
mulate a x^ fitting procedure. Denote by bira{a) the 
spherical harmonics decomposition of a model predic- 
tion for the distribution of gamma-rays depending on a 
set a of parameters describing the model. The best-fit 
parameters can be found by minimizing 



X 



'(«)= E cr,l,m'«^-bUam 



Im, I'm' 



ai'm'-bi'm'ia)), 

(A5) 
where the star denotes complex conjugate and C^^ ;,^, 
is the inverse of the covariance matrix 

/ .^lm.l'm'^l'm',l"m" =6wS,n,n"- (A6) 

I'm' 

Appendix B: Monte Carlo test 

In this appendix we check our method for separating 
a DM fraction by generating random distributions of 
photons. For the test we use an isotropic distribution 
plus a distribution coming from DM annihilation in an 
NFW profile (Equation dH])). 

The model parameters are the same as in Section 
nil Al the window is \b\ > 10°; we use Y/m harmonics 
with / = 0, . . . , 15 and m == 0; the HEALPix parameter 
is nsidc = 32. 
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FIG. 9: Random test of the spherical harmonics decom- 
position algorithm for separating a DM NFW annihilation 
signal from an isotropic background. The input DM frac- 
tion is 0.04. Plotted fractions represent the result of fitting 
in spherical harmonics space. An approximation of these 
fractions by a Gaussian distribution is the "Test results" 
curve. "Test expectation" is a Gaussian around 0.04 with 
the scatter derived from Equation ((Tjl. 



The average total number of photons inside the win- 
dow is lO'^ with an average fraction of photons coming 
from DM annihilation 0.04. In every pixel p we put a 
random number of photons Up according to the Pois- 
son distribution with an average equal to the combined 
density at that pixel ^p — fip°^'^ + ^p ^ . 

In the test we generated iVtests = 100 realizations of 
the photon map. For every realization i, we find the 
best estimate of the DM fraction qf^. The average 
among the realization and the standard deviation are 
(3.98 ± 0.67) X 10~^. The corresponding distribution is 
presented in Figure [S] as "Test results." 

In real applications, there is usually only one real- 
ization of the data available, i.e., we need a way to 
estimate the uncertainty of the result based only on 
one realization. This uncertainty can be estimated 
from the curvature of x^ near the minimum in the di- 
rection of the parameter (Equation (O). The uncer- 
tainty derived from x^ (averaged over the realizations) 
is (7 = 0.68 X 10~^. The corresponding Gaussian dis- 
tribution (4.00 ± 0.68) X 10^2 is plotted in FigureHas 
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"Test expectation." 

The expected deviation of the mean is cr/V-^tcsts ~ 
0.07. Thus the actual deviation of the mean DM frac- 
tion from the expected value is less than one sigma. 
Also the difference of the actual standard deviation and 
the expected one is less than one sigma. We conclude 
that, given a particular distribution of photons, the 
best estimate of the variance given by Equation ([7]) 
is an adequate representation of the actual variance 
among the realizations of the photon distribution. 

The x^ = 13.9±5.2 is calculated from Equation ([M]) . 
The number of degrees of freedom is fourteen: there are 
sixteen data points corresponding to YJo harmonics for 
I — 0, . . . , 15 and two varying parameters correspond- 
ing to the normalization of the two templates: isotropic 
and DM. 

We find that the spherical harmonics decomposition 
is a statistically unbiased fitting method with a viable 
estimation of statistical uncertainty. 



Appendix C: Point sources 

In this appendix we study the dependence of spher- 
ical harmonics on the contribution from point sources. 
We find that in the presence of point sources the vari- 
ance of a;m's increases. 

We will assume some pixelation of the sphere with 
pixels of equal area. Denote by Up the number of pho- 
tons in a pixel p and define the spherical harmonics 
coefficients 



IS 



air. 




(CI) 



where the sum is over the pixels, 7p is the center of 
pixel p, and N-^ is the total number of photons. The 
normalization of a/m's in this appendix is different from 
the normalization everywhere else in the paper. We 
show below that with this normalization the expected 
variance of spherical harmonics in the case of Poisson 
noise is equal to one. 

Analogously to the derivation of the covariance ma- 
trix in Equation (|A3p . we find that the variance of aim's 



A-K 



Var(a,„) = j^Y. >^i™(7p)i1m(7p)Var(n,) (C2) 



If we assume that the diffuse emission and the point 
sources were distributed isotropically, then (a;,„) = 

and (a^„) = (a^n') f*^^ ^^^ ^' "*' "^'^ ^^ ^^^^ case, one 
can relate the variance of spherical harmonics to the 
expectation value of the angular power spectrum 



Var(a;„ 
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1 ' 

— Y. Var(az„) = (G). (C3) 

m— — / 



In the following we will show that this is a good ap- 
proximation for large ^, whereas small (. harmonics are 
dominated by the non-isotropic Galactic emission with 
{aim) ^ 0. 

In the case of the Poisson statistics, the best es- 
timate for the variance of the number of photons is 
Var (rip) = Up. Taking into account that for any point 
on the sphere 



21 



^—^T.YL{i)yui) - 1^ 



(C4) 



we find from Equations (jC2l - IC4I) that in the case of 
the Poisson noise 

Var(a,™) ^ j^Y. Var(np) = 1 (C5) 



Now suppose that there are some point sources. De- 
note by Xm the expected number of m-photon sources 
inside a pixel. In this case the statistics of photons in 
pixels across the sky is not Poisson. Instead, we will 
assume the Poisson statistics of the point sources with 
the average values Xm- In particular the variance of 
the number of photon sources is 



Var(x,nj — Xr, 



(C6) 



The variance of the number of photons in a pixel is the 
sum of the variances of the photon sources times the 
number of photons from every source squared 



Var(np) = y ^ m^Var(a;„i) = y 



TTl XjYi 



(C7) 
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Let us introduce the following parameter 

_ _ Em "^^^m 



(C8) 



where J2m ''^^■m = (rip) is the expected number of pho- 
tons in a pixel. If there is a significant contribution of 
multi-photon point sources to gamma-ray data, then 
'Tiav S> 1, while TOav — 1 for truly diffuse emission. In 
analogy with Equation (|C5p . we find 



Var(o;m) = -7^ X! Var(np) 



Wa 



(C9) 



To summarize, for isotropic distribution of photons 
and for the angular scales smaller than the detector 
PSF, we expect 



iCi] 



(CIO) 



This limit should be saturated for sufficiently large £. 
In the presence of point sources, the variance is TOav 
times larger than the variance in the Poisson statistics 
case. Consequently, for isotropic distribution of point 
sources (or when the angular scale corresponding to i 
is much smaller than the scale of the distribution) , we 
expect 



iCi] 



rria 



(Cll) 



We expect this behavior for intermediate values of £. 
At small £, the Cj's are dominated by a large-scale 
distribution of gamma-rays. 

In Figure 10 we compare the angular power spectra 
before and after masking the gamma-ray point sources 
for several characteristic energy bins. We use the same 
Fermi data as described in Section IlII Al and mask the 
Galactic plane within |6| < 10°. For ^ < 20 the an- 
gular power spectrum is dominated by the large-scale 
distribution of photons: these values of £ are of interest 
for fitting templates at large angular scales using the 
spherical harmonics. For intermediate fs the angular 
power spectrum is dominated by the contribution from 
point sources. For large £ the angular power spectrum 
is consistent with the Poisson noise with an exception 
of the highest energy bin (100 - 300 GeV), where the 
signal in the presence of point sources is above the Pois- 
son noise level even for £ ~ 1000 corresponding to an- 
gular scales ~ 0.2°. This is consistent with the .PSF. 



5- 0.1° for gamma-rays with energies above 100 GeV 

In the analysis we use the HEALPix parameter nside 
— 1024 which corresponds to approximately 10^ pix- 
els. The "pixelized" values of C/'s are supposed to 
be smaller than the values computed by a continuous 
integration, C/pix = wfCi cont, where wi is the pixel 
window function. The window function is a decreas- 
ing function which is equal to 1 at £ = 0. For nside 
= 1024 and £ = 1000, wf w 0.956, i.e. the values of 
Ci's in Figure 10 are less than about 5% smaller than 
the real values. For nside = 32 the window function at 
£ — 15 is wi — 0.989. We did not rescale the spherical 
harmonics by the window function in the fitting proce- 
dure. The only effect of the window function is to put 
a little more weight on lower harmonics. 

Let us make one more technical comment about the 
derivation of the plots in Figure 10. In order to relate 
the variance in the spherical harmonics to the value of 
C/'s we need the average (aim) = 0. In the analysis we 
have used the spherical harmonics decomposition in a 
window |6| > 10° (the window is the same as defined in 
section UlIAp . For the spherical harmonics decomposi- 
tion in a window, at least some of the expected aim 's 
are nonzero. In order to make them zero without affect- 
ing the variance we subtract the spherical harmonics of 
an isotropic distribution inside the window. 
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Angular power spectrum, E = 1.0 - 2.0 GeV 



O 10' 




Angular power spectrum, E = 2.0 - 5.0 GeV 




Angular power spectrum, E = 10.0 - 20.0 GeV 



O 10 




Angular power spectrum, E = 100.0 - 300.0 GeV 




FIG. 10: Power spectrum of angular modes for Fermi gamma-ray data inside a window excluding the Galactic plane 
|6| > 10° for some characteristic energy bins. Blue (upper) lines: all gamma-rays are taken into account, including point 
sources. Red (lower) lines: Fermi gamma-ray point sources [lO|] are masked. The normalization is chosen such that (Ci) = 1 
for the Poisson noise (green constant lines). At low i the power spectrum is dominated by the large-scale structure (the 
red and the blue curves are almost identical). At i '^ 100 the Ci's are dominated by the point sources (the blue curve is 
significantly higher than the red curve). For zero PSF, the blue curve would stay for large i" at a constant level, given by 
Equation (Cll), above the Poisson noise level. The suppression of C;'s to the Poisson noise level for higher £ is due to 
non-zero PSF. In the highest energy bin, 100 - 300 GeV, the PSF is less than 0.1° ll| which corresponds to ^ > 2000. As a 
result, the blue curve stays above the red curve in this energy bin. At lower energies the PSF is higher and the suppression 
of the angular power spectrum to the Poisson noise level occurs for smaller ^'s. In the plots, we use the HEALPix parameter 
nside = 1024. 
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